Random hopping fermions on bipartite lattices: Density of states, inverse 
participation ratios, and their correlations in a strong disorder regime 
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We study Anderson localization of non-interacting random hopping fermions on bipartite lattices 
in two dimensions, focusing our attention to strong disorder features of the model. We concentrate 
ourselves on specific models with a linear dispersion in the vicinity of the band center, which can 
be described by a Dirac fermion in the continuum limit. Based on the recent renormalization group 
method developed by Carpentier and Le Doussal for the XY gauge glass model, we calculate the 
density of states, inverse participation ratios, and their spatial correlations. It turns out that their 
behavior is quite different from those expected within naive weak disorder approaches. 



I. INTRODUCTION 

During the last decade, we have acquired much wis- 
dom of Anderson localization [I], |2( ■ Multifractal scal- 
ing dimensions H for example, have been calcu- 
lated exactly for a random Dirac fermion model [1 , 
and a lot of novel universality classes for disordered sys- 
tems have been discovered 0, H U DU EH D3 • Among 
the universality classes found so far, chiral (or sublat- 
tice) classes H, M have been attracting much interest 
IllllQIIllllllllIl, since they are exactly on ran- 
dom critical points. Gade has predicted for these classes 
that the density of states (DOS) diverges at the critical 
points 0- 

Recently, Motrunich et.al. [19| have claimed that such 
DOS as was predicted within a naive weak disorder ap- 
proach is incorrect if one considers strong disorder fea- 
tures of the model. They have predicted an alternative 
expression of the DOS, which shows a bit weaker diver- 
gence. 

On the other hand, Carpentier and Le Doussal [H l2l) 
have proposed a renormalization group (RG) method to 
analyze strong disorder features at low t emp eratures of 
the XY gauge glass model [H El El El EE E3 Intro- 
ducing higher charge fugacities, they have derived gener- 
alized RG equations for the Coulomb gas model, which 
is remarkably intimate relationship with a nonlinear dif- 
fusion equation called KPP equation |28l l29l y(| . Their 
idea has been applied to a Dirac fermion model with a 
random vector field, and properties of a strong disorder 
regime as well as a weak disorder regime are successfully 
described in a unified way [3ll l32| . 

Interestingly, Mudry et. al. [13 have shown that the 
same RG method can apply to the random hopping prob- 
lem: It has turned out that the disorder strength flows to 
a strong disorder regime, which provides the same DOS 
as Motrunich et. al. have predicted. 

Recent numerical calculations |34j for a random hop- 
ping fermion model on a square lattice with it flux sug- 
gest the divergent DOS, and we can now believe that 
the model actually belongs to the chiral orthogonal class. 
Therefore, we would like to confirm as the next step 
whether strong disorder features of the chiral orthogonal 



class reveal themselves. Unfortunately, the DOS may not 
be a good quantity to grasp them, since the divergent be- 
havior itself is technically hard to observe in numerical 
computations. 

In order to clarify the strong disorder features of the 
model furthermore, we calculate, in this paper, some 
alternative quantities, the inverse participation ratios 
(IPR) and their spatial correlations, where such features 
also appear manifestly. 

This paper is organized as follows: In the next section 
ITT1 we introduce lattice models and their field theory in 
the scaling limit. In Sec. IIHI we derive the scaling equa- 
tions taking into account higher powers of energy terms. 
In Sec. IIVI we convert the scaling equations into KPP- 
type equation and discuss its asymptotic solution. By 
the use of them, we then calculate the DOS, IPR, and 
their spatial correlations in Sec. [V] We give concluding 
remarks in Sec. IVII 



II. MODEL 

In this section, we first introduce the following two 
kinds of lattice models with sublattice symmetry as well 
as time-reversal symmetry: One is the model studied by 
Hatsugai et. al. [l4| defined on a square lattice with it 
flux, and the other is the one on a honeycomb lattice [35| . 
Next we derive a Dirac fermion model in the continuum 
limit of the lattice models at the band center. 



A. Lattice models 

We consider the following non-interacting random hop- 
ping fermions on bipartite lattices in two dimensions: 



H 



(i,3) 



c c, 



(2.1) 



where stands for summation over nearest neighbor 
pairs and tij — to,y + Stij is pure and random hopping 
amplitude, both of which can be chosen as real numbers 
due to time-reversal symmetry. The model defined on a 
square lattice with it flux is here described with a gauge 



2 



to t ij = (— Y y to. In the case of the honeycomb lattice, to,ij 
is uniform, just defined as io.ij = to- 

We are interested in the ensemble-average of the cor- 
relation functions 



= I^WM'Mj)! 292 , 



(2.2) 



where |^„) = £\ *„(i)|i) with 



|0) is a normalized 



wave functions with an energy eigenvalue _E n . Especially 
p(ifi) = p{q) defines the moments of the wavefunction, 
the IPR. To calculate the above correlation functions, the 
following Green functions may be useful: 



1 



r(i)=Im(i|- . 

i — «(w — zi/J 

rfa •*»)(» j) = r« i (i)r« 3 (j). 



(2.3) 



Actually we have 

w 9l+92-l r ( 9 l,?2)(j ;J -) 



91+92 



*n(*)| 29l |^(j)| 292 ^-^), (2-4) 



where numerical constant C q is given by C q = Tr(2q — 
3)!!/(2n - 2)!!. Then, the averaged DOS is given by 
p(E) = r(bo). j n following, we will take the scal- 
ing limit of the lattice model near the zero energy and 
calculate the above correlations by the use of the field 
theory. In that case, it is difficult to define normalized 
wave functions which is crucial to the IPR. To avoid this, 
we alternatively define 



p(?i,92)( 



,_ i) = f^£^M j (2 . 5) 

C91+92 p 



with the normalization condition relaxed. Using these, 
we can define generalized IPR as P^ = p(«>°). Note 
that p in the denominator in Eq. i|2.5[l is needed for the 
"normalization" of the generalized IPR, i.e., P^ 1 ' = 1. 
In what follows, we calculate Eq. (|2.5|l using continuum 
theory and the RG method. 



B. Continuum limit 

In the continuum limit around the zero energy of 
the above lattice models, via ai — > x and Ci/*Ja — > 



J2j e tkFj ipj(x) with a lattice constant a, we have the fol- 
lowing action HI 13 13 HI, after the choice of a suitable 
basis 



S = J cPxr/ji [»7 M (0„ + Af,,) + ihh + M 275 ] fr, (2.6) 

where 71 = <7i, 72 = 02, an d 75 = °3 are usua l Pauli ma- 
trices, and i = 1,2 denotes a "flavor" due to the species 



doubling. The probability distributions of the random- 
ness are 

P[A»] cx exp ("2^ / d 2 xAl 



P[Mu] cx exp 



77^ / d 2 xM 2 u 



(2.7) 



The difference between two lattice models introduced 
above is that the initial strength gA and gm is given by 
(JA = gu for the model on a square lattice while gA 7^ <?M 
for the one on a honeycomb lattice. However, it has noth- 
ing to do with the long distance behavior, since gA only 
is renormalized, as we shall see soon. 
The energy term in Eq. 12.3fl are 



s y = j/i J d 2 xy, 



(2. 



where y-y = u> — iE and y = 0+ + O- with 

O- = 1pL2'ipRl - 1pLllpR2, (2.9) 

Here, chiral fermions have been defined by 



4>Ri 
tpLi 



(2.10) 



In the following section, we add powers of this energy 
term to the action by the use of the replica method. They 
play a crucial role in the RG analysis, as we shall see. 



III. RENORMALIZATION GROUP ANALYSIS 

In the previous work within usual weak disorder ap- 
proaches, the replica method 0, E, 0] and supersym- 
metry method [3 have been used to take quenched av- 
erage. The latter method is quite interesting, since it 
enables us to derive exact beta functions for the present 
model. The recent work by Mudry et. al. has also used 
the supersymmetry method, but with higher energy (or 
"fugacity" in the context of the Coulomb gas model by 
Carpentier and Le Doussal) terms being added, anoma- 
lous dimensions of them may be hard to obtain beyond 
one-loop order. 

Therefore, we will use the replica method for simplicity 
and derive the RG equations up to one-loop order for 
consistency. We follow a similar formulation developed 
by Mudry et. al. j33|. 



A. Scaling equations for fugacities 

Let us introduce the replica, ipi — > ipi a with a replica 
number m; a = 1, 2, — , m. Then ensemble-averaged cor- 
relation functions can be calculated as 



r(9l,92)(a; _ y) 

= (y ai y a2 ■ ■ ■ y aqi (x)y bl y b2 ■ ■ ■ y bq2 (»)), (3.1) 
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in the replica limit m — > 0. 

We utilize the Hodge decomposition to take ensemble- 
average over the vector field, 

= e^d^tp + d^4>, (3-2) 
which converts the probability distribution (|2.7() into 



P[A] oc exp 



1 



d 2 x 



(3.3) 



The gauge transformation 



-11^75 + 



(3.4) 



together with M 1 ±iM 2 -> e ±2lv (Mi±iM 2 ), which gives 
rise to no change to the probability distribution P[M M ], 
yields the replicated action 



with 



J d 2 x I C + A + 2g M M + J2 y™ yn ) '( 3 - 5 ) 



C = 7T— (d^f , 
2gA 



O 



M 



JrolRJ, 



(3.6) 



where we have denoted a = ia and (3 = jb for simplicity, 
and defined currents as 



(3.7) 



Note that the field ip has been decoupled and hence omit- 
ted in Eq. (|3~S|) . 

The last term in Eq. 1)3.5(1 is the nth power of 
the energy term, y n = J2{ a } ' ' ' 3>a» ■ Expand- 
ing these with respect to 0± defined in Eq. I|2.9|l . 
it turns out that the most relevant operators are 

E{a} P+ai ■ ■ ■ +an + O_ 01 ■ ■ • 0_ „) = 0% + 0™ = 

y^ri which is actually converted via the gauge transfor- 
mation l|3.4|) into 



= e 2n ^O r , 



-2ri(f)Qr. 



(3.8) 



As we shall see momentarily, it is important to include all 
these terms even if we started initially with y± only, since 



these terms are generated in the process of the renormal- 
ization. To see this, note that the dimension of the ver- 
tex operator e n ^ is —n 2 gA/^ and therefore, the following 
OPE holds; 



^n</>(z) n'</>(0) 



1 



| z |2nn'g^/n 



,(n+n')<t>(P) 



(3.9) 



This equation is actually involved with a "fusion" of the 
energy terms, generating the higher powers of operators 

y n . 

Based on the above OPE, let us now compute the 
OPE between J^mr an d 3^" r - One problem is the cross 
term such as e"^ z) O^(z)e~"'^ 0) C>I!'(0). Its most singu- 
lar term in the OPE yields the exponent — 2nn'gA/ir + 
\n — n'\, and therefore at least when gA > t/2, such 
OPE gives rise to no singularity. Therefore, we will ne- 
glect these terms, and then we reach the following OPE 
for y n 



3E,(*)3Er(0) 



i 



(3.10) 



This leads to the scaling equations 
dy, 



dl 



(i n = (2 - X n )y n - 7T ^ Vn'Vn-n', (3-11) 



where x n is the scaling dimension of the operator y^ r 
given by x n = n—^-n 2 for the present. However, the Om 
yields an anomalous dimension, which is determined in 
the next subsection. The correlation function r^ qi ' q2 '(x — 
y) is dominated by 



r(*.*)(x -y) = (»))• (3.12) 

In Sec. [V] we will determine its long distance behavior 
using the RG equation for these correlation functions. 



B. Scaling equations for random mass and 
vector-field couplings 

In this subsection, we derive the scaling equation for 
gu- The basic OPE for the currents are 



J 



Jr<xp{z)Jrk\{w) 



3 k 



1 



jRa/3{z)lp 



Rk 



47r 2 (z — w) 2 2tt(z — w) 

^Ra, JRa0{z)lpR K 



[S/3 K JRa\(w) - S a \jR K p(w)] 



2tt(z — w) 



2tt(z — w) 



(3.13) 
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and similar for L-movers. Therefore, we have the follow- 
ing OPE for O m E3 



Om(z)O m (0) 



1 



47T 2 |zl 2 



[2O A (0)-AmO M (0)} 



O M {z)O A {0) ~ O A {z)O A {0) ~ 0, 



(3.14) 
(3.15) 



4tt 2 UI 2 



[y n (0)+A^ n (0)]. (3.16) 



Here, another operator (9^ = JruolJlpp is generated. 
Note, however, that it is just an operator appearing 
after ensemble-average directly over instead of the 
gauge transformation l|3.4l) . and hence it is involved in 
the renormalization of g A ■ Note also that the OPE (|3.16() 
contains operators not included in y n . However, we ne- 
glect them in order to close the OPE algebra. The RG 
equations for g A and gu in the replica limit (m — > 0) are 
then given by 

d 9A _ 9 m 

at 7r 



dg 



M 



dl 



= /3. 



M 



0. 



(3.17) 



and the anomalous dimension of y"l is calculated as 



1 



9M \ 
7T / 



9A 2 

— n . 



(3.18) 



If we neglected the fusion of the energy terms, this would 
serve as the scaling dimension of the operator an d 
govern the long distance behavior of the correlation func- 
tion (y n ). 



It is readily seen that this function G obeys the following 
KPP-type equation, 

^diG = D{l)dlG + G{l-G), (4.3) 

where D(l) = is a diffusion constant. The boundary 
condition is G{— oo,0) = 1 and G(oo,0) = 0. The differ- 
ence from the normal KPP equation is that the diffusion 
constant D depends on /. In analyzing the above equa- 
tion, we take an "adiabatic" approximation as Mudry et. 
al. applied to the same problem. Namely, we first ne- 
glect the /-dependence of the diffusion constant D for a 
while, and after solving the KPP equation, we take the l- 
dependence into account. This scheme may hold at least 
in the leading approximation. It should be noted that 
the function G is a generating functional oiy n - Actually, 
expanding with respect to e x leads to 



gMh)t\ 



(4.4) 



n=l 



In the limit I — * oo, the KPP equation asymptotically 
allows the following traveling wave solution G(x, I) ~ 
G (x — m(l)), where m(l) is a velocity (times I) of the 
traveling wave. In order to determine the dimension of 
the operator 3^ r , let us perturb the action with it, set- 
ting the initial coupling y n = y q S n . q in Eq. (|3.5(l . Then 
it turns out that 

2{Dq+l/q)l + 0{l), for q < 1/VD, 

m (l) = ^ Ay/D(l- §lnZ) +0(1), for g = 1/^4.5) 
4^D (l~ flnZ) +0(1), for q>l/y/D, 



IV. KPP EQUATION 



The average value of u can be evaluated as |2£ 



So far we have derived the scaling equations for y n , 
g A , and gu- Those for g A and <?m are actually scaling 
equations in the replica limit. Basically, we also have to 
solve the equations for y n and take the replica limit. To 
this end, we utilize the method developed by Carpentier 
and Le Doussal |20j . 

Define the following distribution function P(l, u) of y n 



Vn{l) 



7m! 
2 



due nu P(u,l) 



i nu\ 

)p- 



(4.1) 



Then it turns out that the resultant equation for P(u, I) is 
free from m, and we need not to worry about the replica 
limit. Moreover, define 



G(x,l) 



1 



exp 



u—x+(l— gM/ir)l 



(u) P ~m{l)-(l- 9 -f)l, (4.6) 
Using this, we define the typical value of y n as 



2/typ,n(0 = :« 



Ap 



(4.7) 



It is readily seen that yt yP ,n satisfies the following "scaling 
equation" , 



dl 



dm(l) 
dl 

^n2/typ,n i 



9M 



ytyp, 



(4.8) 



(4.2) To be explicit, Eq. jPjl yields 



5 



2- (l- 2M.) n - 
4VD(l-^)n-(l-^) 
4VD(l-^)n-(l-^) 



n, 
n, 



for 
for 
for 



71 < 1/V-D, 

n = 1/-/D, 
n > 1/VD. 



(4.9) 



Notice that the first line of the equation reproduces the (l + dl)a, it turns out that (D^r) obeys the RG equations 



naive dynamical exponent of y n , i.e., z n = 2 — x n , where 
x n is given by Eq. (|3.18(l . We denote the exponent thus 
defined as z w . n , which we shall discuss in Sec. IV CI This 
result implies that effective scaling dimensions of y n can 
be modified due to the fusion of the energy terms. 

Let us now take account of the dependence of D on I 
through gA- It should be stressed that even if one takes 
into account the /-dependence of D in deriving Eq. (|4.9|) 
from Eq. I|4.8|l . the following calculations hold. For suf- 
ficiently large I, z n is given by the last of the selection 
rule above, since 1/y/D ~ 1/V? < 1. Namely, we have 



AnVD 
nz(l), 



z(l) 



const. (J^—j 1 2 ; 



(4.10) 



and using this, we define the typical scaling exponent 
x t yP ,n of the operator 3>™ r as 



x t yP ,n (0 



2-nz(l), (Z->oo). (4.11) 



There are two differences between z n and z WjTl : One is 
the Z-dependence, which is essential to the behavior of the 
DOS, and the other is the n-dependence, which plays a 
role in the spatial correlation functions of the IPR. 



V. THE RG EQUATIONS 

In order to calculate the DOS, IPR, and their spatial 
correlations, we utilize the RG equations for the correla- 
tion functions T^ qi < q2 ^(x — y). Unfortunately, the equa- 
tions are too hard to solve, so that we calculate their 
"typical" values based on the scaling equations for typi- 
cal exponents derived in the previous section. 



A. One-point functions 



In this subsection, we first derive the RG equation for 



r(9>°) to calculate the DOS and the IPR. Recall that 
J^mr includes two fermion correlation functions such as 
(V* Ria( x )4> L2a( x )) > which diverges for a finite u>. In order 
to get finite renormalized correlation functions, point- 
splitting the fields is necessary, Ria^) 1 ^ L2ai x + a ))- 
Then (y^ix)) does not depend on x but does on the 
cut-off a. Therefore under the scale transformation a — > 



where the beta functions are defined in Eqs. (|3.11|) and 
(|3.17l) . and the (matrix of) anomalous dimension 7 is 
given by 



m—q 



n=l 



(5.2) 



with x q defined in Eq. l|3.18[l 6, 36]. Since these equa- 
tions are coupled together, it may be difficult to solve 
and take the replica limit. 

One of possible approximations is to neglect the fu- 
sion of the operators J^mr m Eqs. (|3.10() and JSTHJ. This 
approximation corresponds to the conventional weak dis- 
order approaches. The anomalous dimension ^ q becomes 
in this case trivially diagonal and the RG equation is 
readily integrated. The results thus obtained are sum- 
marized in Sec. IV CI However, this scheme is not able to 
grasp essential strong disorder features the present model 
should have. 

An alternative way is to compute the "typical" values 
of the correlation function, since we have been able to 
calculate yt yp . n in the last section. This includes strong 
disorder effects, as it should be. Assuming that these 
typical values yt yP ,n correspond to the typical values of 
the correlation functions, which we denote as (3^mr)typi 
one can simplify the RG equations above for these typical 
values, since the scaling equations for yt yp ,n in Eq. I|4.8|) 
are "diagonalized" , and so is the matrix of anomalous 
dimension j n . The RG equations for the typical values 
should be 

dl ~ &A dg~A + XtyP V {yL)tyP = °" (5 ' 3) 
We are then led to 



(JC-)typ ~ ex P 



dl' 



; (0 



(5.4) 



On the other hand, recall that the energy E is related to 
?/i as y\ = uj — iE. Corresponding to ytyp.i, let us define 
the typical energy yt yp ,i = uit yp — iEt yp . From Eq. (|4.8J) . 
it follows that 



A 



typ 



typ 



exp 



exp 



dl'z{l') 



const 



(9_M\ 



(5.5) 
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where A typ is a renormalized energy and const, denotes Therefore, we reach the following IPR: 
a nonuniversal positive constant. Hence, as a function 



of Btypj we can rewrite the typical DOS ptyp = (3^mr) 



t-yp 



into 



Ptyp(23typ) 



A 



typ 



exp 



typ 



In- 



A 



typ 



typ 



(5.6) 



P^JE) ~ const. 



(5.10) 



This implies that the multifractal scaling exponent r(q), 
defined by P^ ~ L~ T ( q \ is given by 



with 



const. I 

,9m J 



This has been derived for the first time by Motrunich et. 
al. [Tgj with the help of a str ong coupling expansion, and 
rederived by Mudry et. al. [M] by using the supersym- 
metry method. What is responsible for this value of k 
is the velocity selection rule (I4.5|l or (|4.9() and resultant 
dynamical exponents l|4.10|l . Recent numerical calcula- 
tions [3^| support, to be sure, the divergent behavior of 
the DOS, but the nonuniversal constant in Eq. I|5.7() may 
be too small for us to observe the exponent k in the nu- 
merical calculations, unfortunately. 

To overcome the difficulty, we calculate the IPR for 
general q, which could make the strong disorder features 
manifest. For general q, we have 



w typ Omr)typ 



Ptyp(0. 



T(q) = 0, 



(5.11) 



(5.7) for the typical values. This result is also attributed to 
the velocity selection rule (|4.9[1 and therefore the strong 
disorder features of the present model. For references, 
we summarize in Sec. IV CI the results of a naive weak 
disorder approach. 



B. Two-point functions 

Constant IPR obtained in the last subsection usually 
implies a localized wavefunction. However, it is believed 
in general that the zero-energy wavefunction of chiral 
models is extended. To clarify this point, we calculate 
the spatial correlations of the IPR in this subsection. The 
correlation function T^^\x - y) = (y^(x)y^r{y)) 
(5.8) obeys the following OPE 



which follows from 



(q - l)z + Xtyp.q 



(q-l)z+(2-qz) 
2- z 

£typ,l- 



(y&(x)y%Ay))~C qiq2 (x-y){y^(y)). (5.12) 
This may be valid for a <C \x — y\ <C L, where L is 
a system size, and we will determine the r = \x — in- 
dependence of the coefficient C qiq2 . It is easily verified 
(5.9) that the coefficient C qiq2 obeys the RG equations, 



d d d 

r - Pn-R PA-r, — 

dr dy n dg A 



Cgi<22( r ) C(Ji<72 ( r )'lq 1 +q2 



(^ 1 I +92 > + (7 gi +7 92 )(n i r n 2 r )=0, 



(5.13) 



where 7 is defined in Eq. I|5.2() . These coupled set of 
equations become tractable for the typical values: 

d d \ 

t-q^ - Paq^ + Xt yp><n>i2 J c typ,nq2( r ) = °. ( 5 - 14 ) 

where Xt yPtqitq2 = x typ , g i +^typ, 92 - Xtyp, gi+(?2 .I^r suffi- 
ciently large I, it follows from Eqs. H4.10[) and 14.1111 that 



X, 



typ, gi 1 qi 



2. Thus we obtain 



p(*?l,<?2)/ 

^typ \ x 



y) 



const. 
\x - y\' 



(5.15) 



where we have used that previous result l|5.10fl . Remark- 
ably, the result shows that the exponent is a constant, 2, 
independent of q\ and q 2 as well as the disorder strength. 
Since the correlation functions show power-law behavior, 
the wavefunction should not be localized. 



C. Weak disorder approximation 



So far we have obtained several formulas taking into 
account strong disorder effects. As has been discussed 
by Motrunich et. al. as well as Mudry et. al., the ex- 
ponent k in the DOS is k = 2/3, which is different from 
K = 1/2 obtained previously by Gade with the help of 
a naive weak disorder approach 0- In order to clarify 
the strong disorder features of the IPR and their spa- 
tial correlations obtained in the previous subsection, we 
calculate them within the conventional weak disorder ap- 
proach, for reference. 

Without considering the fusion process of higher en- 
ergy terms, the dynamical exponents of y n is given just 



7 



by the first line of Eq. I|4.9(l . Namely, 



2Dn 2 = n 2 z w (l), 
~> const. (^M-^j I. 



(5.16) 



Here, the label w means exponents within a weak disorder 
approximation. Note the difference of the /-dependence 
between Eqs. (|4.10|l and (|5.16() : This reflects the expo- 
nent k in the DOS. Neglecting the fusion of the operators 
y™ r , we find that the RG equation for (3^mr) is diagonal, 
given by Eq. (|5.3[) but with the naive scaling dimension 
x q ~ 2 — q 2 z w (l) in Eq. (|3.18() . Therefore, together with 



A 

E ^ 6XP 



exp 



dl'z w {l') 
cons t. (iM.) f 



\ 7T / 

we have the same DOS (15.611 but with 

_ 1 

Kw - 2" 



(5.17) 



(5.18) 



Next, let us consider the correlation functions for gen- 
eral q. The solution of the RG equations yields 

~ exp j- J dl' [(q - l)z w - x q + xi]\ 



= exp 



dl' z^ 



(5.19) 



from which, together with Eq. (|5.17() . we reach 



-9(9-1) 



(5.20) 



Negative exponent in the above equation suggests that 
the present weak disorder approach is incorrect and the 
result in the last subsection (|5.10(l is more reasonable. 
Since we are dealing with the same model in the same 



formalism, this result also implies that the exponent k in 
the DOS should be k = 2/3 rather than 1/2. 
In a similar way, we reach, 



exp 



-cqiq 2 (ln|x - y\Y 



(5.21) 



where c is given by the same disorder-dependent constant 
as in Eq. I|5.7|l . 



VI. CONCLUDING REMARKS 

In this paper we have reconsidered random hopping 
fermion models, based on the recent developments of the 
RG method for a strong disorder regime. We have used 
the replica method and calculated the DOS, IPR, and 
their spatial correlations. We have been able to repro- 
duce the DOS predicted by Motrunich et. al. as well as 
Mudry et. al. Moreover, we have calculated the IPR and 
their spatial correlations, whose behavior near the zero 
energy is also different from the ones predicted by using 
conventional weak disorder approach. 

In the present calculations, including higher powers of 
the energy perturbation (or higher fugacity terms in the 
context of Coulomb gas model) plays a crucial role, which 
has been invented for the first time by Carpentier and Lc 
Doussal and applied to the present model by Mudry et. 
al. 

Although recent numerical calculations have supported 
the divergent DOS toward the zero energy, it may be 
difficult to distinguish the exponent k. We hope that 
IPR and their correlations calculated in the present paper 
may be useful to establish the strong disorder features of 
the random hopping models. 
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